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A CUT DISCONTINUOUS GALERKIN METHOD FOR THE 
LAPLACE-BELTRAMI OPERATOR 


ERIK BURMAN, PETER HANSBO, MATS G. LARSON, AND ANDRE MASSING 


Abstract. We develop a discontinuous cut finite element method (CutFEM) for the Laplace- 
Beltrami operator on a hypersurface embedded in R d . The method is constructed by using 
a discontinuous piecewise linear finite element space defined on a background mesh in 
The surface is approximated by a continuous piecewise linear surface that cuts through the 
background mesh in an arbitrary fashion. Then a discontinuous Galerkin method is formulated 
on the discrete surface and in order to obtain coercivity, certain stabilization terms are added on 
the faces between neighboring elements that provide control of the discontinuity as well as the 
jump in the gradient. We derive optimal a priori error and condition number estimates which 
are independent of the positioning of the surface in the background mesh. Finally, we present 
numerical examples confirming our theoretical results. 
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1. Introduction 

1.1. Background. Recently there has been a rapid development of so-called cut finite element 
methods (CutFEM), which provide a technology for discretization of both the geometry of the 
computational domain and the partial differential equations (PDE) which we seek to solve. The 
basic idea in CutFEM is to represent the geometry of the domain on a fixed background mesh 
which is also used to construct the finite element space. The geometric representation of surfaces 
consists of elementwise smooth surfaces that are allowed to cut through the background mesh in 
an arbitrary fashion. The active background mesh then consists of all elements that are cut by the 
domain and the finite element space used in the computation is the restriction to the active mesh. 
The variational formulation of the PDE is defined on the cut elements. Boundary and interface 
conditions are enforced weakly. In order to obtain a stable method, independent of the position 
of the geometry in the background mesh, and to handle the cut elements in the analysis, certain 
stabilization terms are added that provide control of the local variation of the discrete functions. 
Further stabilization may be necessary to control for instance convection or to establish an inf-sup 
condition. CutFEM can be used to handle problems in the bulk domain, on surfaces, and coupled 
bulk-surface problems. 

1.2. Earlier Work. While the development of standard finite element methods on triangulated 
surfaces for the numerical approximation of surface PDEs was already initiated by the seminal 
paper by Dziuk m, CutFEM-type methods for surface PDEs have been introduced only recently. 
Olshanskii et al. [32 proposed the first discretization of the Laplace-Beltrami operator based 
on restricting continuous piecewise linear finite element functions from the ambient space to the 
surface. The matrix properties of the resulting system were investigated in Olshanskii and Reusken 
m showing that diagonal preconditioning can cures the discrete system from being severely 
ill-conditioned. As an alternative, stabilization techniques based on face stabilization and full 
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gradient stabilization were introduced and analyzed by Burman et al. [B] and Deckelnick et al. 
m, respectively. Face stabilization techniques were also used in Burman et al. !?] to develop a 
cut finite element method for stationary convection problems on surface, while Olshanskii et al. 
[33) employed streamline diffusion techniques. For evolving surface problems using CutFEM-type 
techniques, we refer to Olshanskii et al. [34] and Hansbo et al. [25]. Higher order methods can 
be found in Grande and Reuskcn m and adaptive methods including a posteriori error estimates 
were presented by Demlow and Olshanskii tig and Chernyshenko and Olshanskii 0. 

For bulk problems we mention the following contributions: Interface problems were considered 
in Hansbo and Hansbo [22] with similar techniques being used in Hansbo et al. 231 and Massing 
et al. [29] to develop overlapping mesh methods. Face-based so-called ghost penalties were then 
employed to solve the Poisson boundary problem [2], the Stokes boundary problem HIM! and 
Stokes interface problems [MUST]. Alternative CutFEM schemes for the Stokes interface problem 
can be found in Grofi and Reusken m where the pressure space were enriched in the vicinity of 
the interface, and in Heimann et al. 2^, Sollie et al. [3Sj which are based on unfitted discontinuous 
Galerkin methods using cell-merging techniques problems to obtain stable and well-conditioned 
numerical schemes. Higher order discontinuous Galerkin with extended element stabilization for 
an elliptic problem were investigated by Johansson and Larson m- For coupled surface-bulk 
problems, see Burman et al. [8], Gross et al. [2D] , and Massing et al. [25] for implementation issues. 
We refer to the overview article Burman et al. [3] on cut finite element methods and the references 
therein. Finally, we want to mention Dziuk and Elliott [14] for a recent and comprehensive 
overview on methods for surface PDEs. 

1.3. New Contributions. In this paper we develop a stable cut discontinuous Galerkin method 
for the Laplace-Beltrami operator on a hypersurface in R d . Discontinuous Galerkin methods gener¬ 
ally provide excellent conservation properties and advantageous flexibility in terms of mesh refine¬ 
ment and local space enrichment and additionally exhibit good stability properties for convection- 
diffusion problems and problems involving discontinuous coefficients. Extending the results from 
Burman et al. [6], this paper is to the best of our knowledge the first step in developing CutFEM- 
type discontinuous Galerkin method for surface problems. For triangulated surfaces a discontin¬ 
uous Galerkin method was proposed and analyzed in Dedner et al. m The analysis of the cut 
discontinuous Galerkin method on surfaces however poses additional challenges and thus requires 
different technical tools. 

We consider discontinuous piecewise linear elements on a background mesh consisting of sim- 
plices in R d . The discrete surface is continuous and is on each element given by a plane that cuts 
through the element. Such a surface may, e.g., conveniently be constructed by taking the zero level 
set of a continuous piecewise linear approximation of the distance function. We assume that the 
discrete surface converges to the exact surface in such a way that the error in position and normal 
is 0(h 2 ) and O(h), respectively. On the piecewise linear surface we formulate a discontinuous 
Galerkin method. 

In order to prove basic stability results the formulation must be stabilized and we use consistent 
terms that provide control of the variation of the finite element functions over the faces in the 
background mesh. This is achieved by controlling the jump in the function as well as the gradient 
across interior faces in the active mesh. Essentially, the stabilization terms enable control of 
functions on an element in terms of neighboring elements, which allows us to handle the presence 
of small and deteriorated surface elements. 

Two technical estimates are critical to our analysis. First, an inverse estimate of the co-normal 
fluxes at edges in terms of the tangent gradient and the stabilization terms and secondly, a discrete 
Poincare estimate for discontinuous piecewise linear functions that provides control of the L 2 norm 
on the active mesh in terms of the tangent gradient and the stabilization terms. These results 
extend the corresponding results in Burman et al. [B] to discontinuous piecewise polynomials and 
also provide certain improvements in the details of the proof. Based on the stability results we 
derive a priori estimates of the error in the energy and L 2 norm that takes both the approximation 
of the solution and the aproximation of the geometry into account. Furthermore, again using the 
stabilization, we prove an optimal upper bound for the condition number. We emphasize that 
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all our results are completely independent of the relative position of the discrete surface in the 
background mesh and only information available from the discrete surface is used in the definition 
of the method. 

1.4. Outline. In Section [2] we present the model problem and formulate the numerical method, 
while in Section [3j estimates related to the error resulting from the approximation of the hyper¬ 
surface are collected. Section [4] is devoted to the construction of an interpolation operator and the 
formulation of the necessary interpolation estimates. Afterwards, we develop stability estimates 
for the proposed numerical method in Section [5j including certain Poincare estimates for finite 
element spaces on thin domains. We also prove novel inverse estimates for the co-normal flux 
accounting for irregular surface geometry discretizations which are typical in CutFEM methods. 
The core of Section [6] derives a priori error estimates based on a Strang Lemma approach, followed 
by providing bounds for condition number in Section [7] Finally, the numerical studies presented 
Section [8] serves to corroborate our theoretical findings and to illustrate the importance of the 
employed stabilization techniques. 


2. Model Problem and Finite Element Method 


2.1. Preliminaries. In what follows, T denotes a compact and oriented C^-hypersurface, k ^ 2 
without boundary which is embedded in and equipped with a normal field n : T —> R. d of 
class C k ~ 1 . We let p £ C k (U$ 0 (T)) be the signed distance function induced by the normal field n 
and uniquely defined on the tubular neighborhood Ug 0 (T) = {x £ : dist(a;,r) < d 0 } for some 

<5q > 0, see Gilbarg and Trudinger m- For such a tubular neighborhood, p{x) denotes the closest 
point projection given by 

p{x) = x — p(x)n(p(x)) (2-1) 

which maps x £ Ug 0 (r) to the unique point p{x) £ F such that | p(x) — x\ = dist(:r, T). The closest 
point projection allows to extend any function on T to its tubular neighborhood t/, 5 0 (r) using the 
pull back 

u e (x) = uop(x) (2.2) 

In particular, we can smoothly extend the normal field nr to the tubular neighborhood t/< 5 0 (F). 
On the other hand, for any subset T C Us 0 {T) such that p : T —> T is bijective, a function w on T 
can be lifted to F by the push forward 

(w l (x)) e = w l op = w on T (2.3) 


A function u : P —> R is of class C*(r), l ^ k if there exists an extension u £ C l (U ) with it|r = u 
for some d-dimensional neighborhood U of T. Then the tangent gradient Vr on T is defined by 


Vru = PrVu 


(2.4) 


with V the gradient and Pr = Pr(x) the orthogonal projection of onto the tangent plane 
of T at x £ r given by 

Pr = / — nr <S> Rr (2-5) 


where I is the identity matrix. It can easily been shown that the definition (2.4) is independent of 
the extension u. We let ||u>||r = {w,w )r denote the L 2 (T) norm on T and introduce the Sobolev 
H m (T) space as the subset of L 2 functions for which the norm 


IMIrn,r = E \\ D r k M\r, m = 0,1,2 (2.6) 

fc=o 

is defined. Here, the L 2 norm for a matrix is based on the pointwise Frobenius norm, D^'°w = w 
and the derivatives Dp 1 = PpVw, Dp’ 2 w = Pr(V ® Vu>)Pr are taken in a weak sense. Finally, 
for any function space V defined on T, we denote the space consisting of extended functions by 
V e and correspondingly, we use the notation V 1 to refer to the lift of a function space V defined 

on r. 
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2.2. The Continuous Problem. We consider the following problem: find u : T —» R such that 

-A r u = / on T (2.7) 

where Ar is the Laplace-Beltrami operator on T defined by 

A r = V r • V r (2.8) 

and / G L 2 (T) satisfies f r f = 0. The corresponding weak statement takes the form: find u G 
i7 1 (F)/R. such that 

a(u,v)=l(v) \/v € H 1 ( r)/K (2-9) 

where 

a(u,v) = (V r u, V r v)r, l(v) = (f,v) r (2-10) 

and (v,w)r = f r vw is the L 2 inner product. It follows from the Lax-Milgram lemma that this 
problem has a unique solution. For smooth surfaces we also have the elliptic regularity estimate 

Nkr<||/||r (2.11) 

Here and throughout the paper we employ the notation < to denote less or equal up to a 
positive constant that is always independent of the mesh size. The binary relations > and ~ are 
defined analogously. 

2.3. The Cut Finite Element Space. Let Th be a quasi uniform mesh, with mesh parameter 

0 < h < ho, into shape regular simplices of a open and bounded domain fl in containing 

Ug 0 (r). On Th- let ph be an continuous, piecewise linear approximation of the distance function p 
and define the discrete surface as the zero level set of ph', that is 

ri, = {ie!l: Ph{x) = 0} (2.12) 

We note that Tf,, is a polygon with flat faces and we let rih be the piecewise constant exterior 
unit normal to We assume that: 

• T/j C U So (r) and that the closest point mapping p : F^ —► T is a bijection for 0 < h < h 0 . 

• The following estimates hold 

IMU°°(r h ) ^ h 2 , \\n e - n h \\ Lao( r h) < h (2-13) 

These properties are, for instance, satisfied if ph is the Lagrange interpolant of p. For the mesh 
Th, we define active background mesh 

T h = {Tef h -.Tn r h ^0} (2.14) 

and its set of interior faces by 

Th = {F = T + (~l T~ : T + ,T~ £ Th} (2.15) 

The face normals rip and Up are then given by the unit normal vectors which are perpendicular on 
F and are pointing exterior to T + and T~ , respectively. The corresponding collection of geometric 
entities for the surface approximation Th are denoted by 

K h = {K = T h nT:T€%} (2.16) 

£ h = {E = K+ n K~ : K + , K~ G T h } (2.17) 

To each interior edge E we associate the normals rip given by the unique unit vector which is 
coplanar to the surface element K ± , perpendicular to E and points outwards with respect to K ± . 
Note that while the two normals rip only differ by a sign, the normals rip do lie in genuinely 
different planes. The various set of geometric entities are illustrated in Figure [I] We observe that 
the active background mesh Th gives raise to a discrete or approximate /i-tubular neighborhood 
of T h , which we denote by 

N h = U T £T h T (2.18) 

Note that for all elements T G Th there is a neighbor T' G Th such that T and T' share a face. 
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Figure 1. Domain set-up 


Finally, let 

X4=0 A CO (2.19) 

TeT h 

be the space of discontinuous piecewise linear polynomials defined on Th and define the subspace 

V h ,o = {v G V h : Xr h (v) = 0} (2.20) 

consisting of those v G 14 with zero average Ar h (u) = f r v. Note that 14,o can be considered as 
the discrete version of H 1 ( r)/R. 

2.4. Discontinuous Galerkin Method. Before we formulate the cut discontinuous Galerkin 
method for the Laplace-Bcltrami problem (2.71, we introduce the notation of average and jump 
first. As the co-normal vectors are generally not collinear, the average flux for a vector-valued, 
piecewise continuous function cr on A4 is defined by 

{^)\e = (Te ■ n % - (Te ■ n E ( 2 - 21 ) 

while for any, possibly vector-valued, piecewise continuous function w on A4, the jump across an 
interior edge E G O is defined by 

Mb = We ~ w e (2.22) 

with w(x) ± = lim f _ >0 + w ( x ~ trig). Similarly, the jump across an interior face F € Eh is given by 


Mb = wj- w F . 

Next, we define the bilinear form A^(-, •) by 

A h (v,w) = a h (v, w) +jh{v,w) \/v,w G 14 

with 


(2.23) 

(2.24) 


a h (v,w) 


jh(v,w ) 


^ (yr h v,S7r h w) K 

K&K h 

- ((Vu),M)^- y (M>( Vw ))^ 

E£_Eh E^Eh 

+ X] p E h- l ([v\, [w])e 

E€Eh 

P F h~ 2 ([v], M)f +l{n F ■ [S7v],n F • [Viu])f 


(2.25) 

(2.26) 
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where are positive parameters. The right hand side is given by linear form 

k{v) = (f e ,v) r h 

Then the cut discontinuous Galerkin finite element method for the Laplace-Beltrami problem ( |2.7| 
takes the form: find Uh G 14,o such that 


(2.27) 


A h (u h ,v) = l h (v) \/v G T4,o 


(2.28) 


We conclude this section by introducing suitable norms for the forthcoming stability and error 
analysis. Here and throughout this work, any norm || • || -p h involving a collection of geometric 
entities 74 is defined in the usual way, i.e., by summing over all entities: || • ||p h = Y^p^v, II ' lip- 
For the discrete energy norm associated to the weak formulation (2.28) we chose 

IIMIIfc = \\^r h v\\ 2 Kh + \\h- 1/2 [v]\\l h + ||/i _1 Hllk + II n F • [Vu]||^ 

and for the function space ( H 2 (T)) e + 14, we define 

-\\h 1/2 (Vr h v)\\ 2 £h 


\*,h 




(2.29) 


(2.30) 


For convenience, the part of the energy norm associated with the face terms is summarized by the 
notation 

IIMIIk = II^MIIk. + ll™P • l^ v ]\\p h (2.31) 

Although it is not obvious from the definition, we will later prove in that the discrete energy 


norm (2.29) indeed defines a norm on 14,o, see Corollary 5.27 


Remark 2.1. We point that thanks to inverse estimate (4.8), the edge-based, penalty term appear¬ 
ing in (2.25) can be controlled by its face-based counterpart in the stabilization form (2.26). More 

lk“ 1/2 Mlk < II^MIU 


precisely it holds 


(2.32) 


and consequently, one might simplify formulation (2.25) by setting /3e = 0 at the expanse of 
possibly requiring a larger value for ftp. 

2.5. Lifted Discontinuous Galerkin Form. The mismatch of the smooth surface T and its 
discrete counterpart T h gives raise to a geometric variational crime in the proposed discretization 
scheme which must be accounted for in the forthcoming a priori analysis. An important instrument 
in assessing the error introduced by this variational crime will be the following lifted version of 
the bilinear form (2.25): 


l (v,w)= Y (V r u,Vr w) K i 


(2.33) 


K l eic{ 


- Y (( v ^).M)s‘- Y (M AVw)) e i 

E l es‘ h E‘e£ l h 

+ Y 

E‘e£L 


[w]) E l 


(2.34) 


where, referring to the notation in Section 


2.4 


the average (•) is taken with respect to ,( x ) £ 


T x K l A- that is the unique co-normal vector which points outwards K 1 ^ and is orthogonal to both 
the surface nr(x) and the tangential space T x E l,± . With the definition of the lifted discontinuous 
Galerkin form, it is clear that the solution u of problem (2.7) satisfies the following weak problem: 


a l h {u,v) = l(v) Vn G Vjj 

For v G H 2 (T) + V ] h , the natural energy norm associated with this weak form is given by 


kill* = llvHIjk + Ik 1/2 M 


\\h 1/2 (V r v)f £l 


(2.35) 


(2.36) 
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3. Domain Perturbation Related Estimates 

The aim of this section is to collect the appropriate geometric estimates which are necessary to 
quantify the error introduced by the geometric variational crime, i.e., the piecewise linear approxi¬ 
mation of the smooth surface T. While developing face and edge-related estimates simultaneously, 
we essentially follow the presentations in Dziuk Dziuk and Elliott nn. Olshanskii et al. [3^1 
and include without detailed proofs those results which are well known. Here and throughout the 
remaining work, we write (•, -) R d and || • || R d to denote the standard scalar product and its associate 
norm in 


3.1. Gradient of Lifted and Extended Functions. Derivatives of extend and lifted functions 
can easily be computed by the chain rule, once the derivative of the closest point projection is 
known. Using the fact that Vp = nr, the derivative of the closest point projection computes to 

Dp = I — nr <8> nr — pV (g> Vp (3-1) 


where 77 = Vi8)Vp = Vnr denotes the Hessian of the signed distance function. The symmetry 
of the Hessian and the projection P r together with simple fact that 0 = V||n r ||^ d = 2Vnpn r 
implies that Hnr = 0 and (nr <8> nr)H = 0 and therefore that HPr = H = PrH. This allows to 
rewrite 13.11 as 

Dp = P r (I — pH) =P T -pH (3.2) 

and thus 


Dv e = D{v op) = DvDp = DvP r (I - pH) (3.3) 

Exploiting the self-adjointness of Pr, Pr h , and H 1 we have for any vector a £ R. d 
(^r h v e ,a) K d = (Vu e ,Pr h a) R d 

= Dv e Pr h a = DvP r (I - pU)Pr h a = (Vn,P r (J - pP)Pr h a (3.4) 

= (Pr h (I-pH)P r Vv,a) R d (3.5) 

that is 


Vr h u e = P T V r u (3.6) 

where we introduced the invertible linear mapping 

B = P r (7 - pH)P rh : T X {K) -> T p(x) (T) (3.7) 

which maps the tangential space of K at x to the tangential space of T at p(x). Setting v = w l 
and using the identity ( w l ) e = w, we immediately get that 

\7 r w l = B~ T S7 Th w (3.8) 


for any elementwise differentiable function w on T^, lifted to T. From its definition (3.7), it is easy 
to derive norm bounds for operators involving B once H can be controlled. We recall from HU 
Lemma 14.7] that for x £ Us 0 (T), the Hessian H admits a representation 


d 

n ( x ) = 

i= 1 


1 + p{x)nf 


aP ® a,- 


(3.9) 


where iq are the principal curvatures with corresponding principal curvature vectors cq. Thus 


l|ft|U~(L7, o( r)) < 1 (3-10) 

for <5 0 > 0 small enough and as an almost immediate consequence, we have the following estimates 
summarized in 


Lemma 3.1. It holds 
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Proof. The proof is standard, see Dziuk and Elliott m, and only sketched here for complete¬ 
ness. The first two bounds follow directly from (3.7) and (3.10) together with the assumption 
||p||l“(i\) < h 2 . Similarly, it follows that Pr — BB 1 = Pr — PyPr h Py + 0{h 2 ). An easy calcu¬ 
lation now shows that Py — PyPy h Py = Py(Py — Py h ) 2 Py from which the desired bound follows 
by observing that Pr — Py h = (n — rih) 0 n + 0 (n — rih) and thus ||(Pr — Pr fe ) 2 


\n - n h 


'(Th) 


<h 2 


li”(r h 


3.2. Change of the Integration Domain. Next, we derive estimates for the change of the 
Riemannian measure when surface and edge integrals are lifted from the discrete surface to the 
continuous surface and vice versa. For this purpose, we define the quotient of measures |P|d-i 
and |P|d -2 by 

d K l = |P| d _i d K, d E l = |P| d _ 2 d E (3.12) 


where d I\ and d E denotes the measure on K £ K-h and E £ £h, respectively, while the corre¬ 
sponding measures on the lifted entities K l £ K, l h and E l £ £' h are denoted by d K l and d E l , 
respectively. Picking an element K and one of its boundary edges E together with its outer co¬ 
normal he, we can assume (after some rigid motion) that I\ C K d_1 x {0} and E C R d_2 x {0,0} 
and that nx = e<i and tie = Sd-i- In this coordinate system, we have dA' = dx\ ... dxd-i and 
d E = dxi... dxd-2- 

Recall that for any smooth parametrized fc-dimensional sub-nranifold M C together with 
a smooth parametrization p : U C satisfying p(U) = M, the integral f M f d M of a 

function f : M M. is given by 

[ fdM= [ f(p(x))Jdet(g i:i (x))dx, (3.13) 

J M JU V 

where gij = (diP,djp) R d i,j = 1 is the first fundamental form of M in local coordinates. 

Considering the cases M = K l , U = K and M = E l and U = E with the closest point projection 
p as parametrization, we see that |P|d_i and |P|d _2 are precisely given by the d — 1 and d — 2 
volumes 


I B\ d = ^det 

\B\d-2 = fdet ({gij))^ 


(3.14) 

(3.15) 


Instead of calculating these volumes exactly, we derive a simple asymptotic representation in h. 
Combining the representation (3.2) of Dp with the fact that the choice of our coordinate system 
implies 


rip = (rip, ei) R d < h, i = 1,..., d — 1, and thus rip ~ 1 + 0(h 2 ) (3.16) 


allows us to derive an asymptotic expression for the coefficients gtj for i,j = 1,..., d — 1: 

9ij = ((Pr - pH)ei, (Py - pB)ei) R d (3.17) 

= (Pye il Pye l ) R d + 0(h 2 ) (3.18) 

= Sij — riyn J r + 0(h 2 ) (3.19) 

~ Sij + 0(h 2 ). (3.20) 


Consequently, ^Jdet(gij) ~ yj 1 + 0(h 2 ) ~ 1 + 0(h 2 ) and thus 

|B|d~ l + 0(h 2 ) 
\B\d-l ~ 1 + Q(h 2 ) 


(3.21) 

(3.22) 


which leads directly to the bounds summarized in the next Lemma. 


Lemma 3.2. The following estimates hold 

111 - \B\d-l\\L°°{K h ) < h~, |||P|d||z,“(/C h ) ~ 1) 

111 - |-B|d-2||.L'»(£ h ) ^ B, |||S| d _ 2 ||Loo(f h )< 1, 


Ill-Bid 1 \\l°°(k h ) ^ 1 
\\\B\- d l \\ L ~ {£h) <1 


(3.23) 

(3.24) 
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We point out that estimates similar to (3.24) were derived in the finite volume context in 
Giesselmann and Muller m and used in Dedner et al. El- 


Next, we note that by combining the estimates for the operator norms (3.11) and for the metric 
distortion factors (3.23) ( [3.24 ), the following norm equivalences can be obtained: 

Lemma 3.3. Let v £ L 2 (Vh) and w £ L 2 (V l h ) for Vh £ {ICh,£h}- Then it holds 


Mv h , 




\\ w \\Vh 


If in addition v £ and w £ the following equivalences are satisfied 


||ViV 


\Vr h v\\ 


HVrHI 


llVi^Xlk, 


(3.25) 


(3.26) 


Before we proceed, let us note that the measure quotients (3.14)- (3.22) can be equivalently 
expressed using the exterior product of the tangential vectors Opp. i = 1,. .., d — 1 and the surface 
normal rip. More precisely, 

|B| (l _ 1 = ||&pA...A3 d _ip|| (3-27) 

\B\ d - 2 = \\d x p A ... d d _ 2 P A n r || (3.28) 

,v d -i £ the outer product v\ A ... A v d -± is the unique vector 


where for given vectors v \,. 
satisfying 

det(vi, 


, v d -i,w) = (v\ A ... A v d -i, w) R d \/w £ 


(3.29) 


With this in mind we can now easily estimate the deviation of the lifted co-normal vector n l E from 
the co-normal vector n E i associated with the lifted edge E l . 

Lemma 3.4. It holds 

II B l n l E — \B\ d _2n E t \\l°°(e 1 ) ^- 2 , (3.30) 

Proof. Expanding the scaled co-normal and lifted discrete co-normal vectors in the standard or¬ 
thonormal basis ei,... ,e d , it suffices to show that the resulting expansion coefficients satisfy the 
estimates 

( B l n l E - \B\ l d _ 2 n E i ,e») K d < h 2 i = l,...,d. 

To compute ( B l n E ,ei)s.d , we simply calculate 


(3.31) 


(■B n E , ei) R d — ( Pt{I ~ pTL)e d - 1 , ej) R d — (Pre d _i, e,) R d -I- 0(h ) 


= S d - lti - nf-X + 0(h 2 ) 

(0(h 2 ), i=l,...d-2 

= <l + 0(/i 2 ), i = d- 1 
[— np _1 rip + 0(h 2 ), i = d 


(3.32) 

(3.33) 

(3.34) 


Now turning to — (\B\ l d _ 2 n E i , ej) R d, we observe that due to (3.28), the fact that — n l E is pointing 
inwards and n E _L d\p, ..., d d _ 2 Pi n r, we have 

-\B\ d _ 2 n E i = dip A • • • A d d - 2 p A n r , (3.35) 


and therefore by definition of the exterior product (|3.29 1 
— {\B\d- 2 n E i,e, 


= det(Dpex ,..., Dpe d - 2 , nr, ef) 

= det(P r ei ,... ,P r e d - 2 ,n r ,ei) + 0(h 2 ) 

= det(ei,..., e d - 2 , nr, ef) + 0(h 2 ) 

0 + 0(h 2 ), i = l,...,d — 2, 

= —np + 0(h 2 ) = — 1 + 0(h 2 ), i = d — 1, 


(3.36) 

(3.37) 

(3.38) 

(3.39) 


d -1 


+ 0(h 2 ) i = d. 


Now adding (|3.34|) and (|3.39|) we arrive at the desired estimate by observing that np 1 — nf 1 n < f = 
np _1 (l — np) = Rp -1 ■ 0(h 2 ). □ 
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The previous lemma roughly states that push-forwarding a properly scaled discrete co-normal 
vector results in a good approximation of the co-normal vector field along the lifted edge. As a 
result we can prove the final lemma of this section. 

Lemma 3.5. For v £ H 1 (T) e + Vh it holds 


< 




(3.40) 


Proof. We only sketch the proof. Recalling the norm definitions (2.29) and (2.36) and the norm 
equivalences from Lemma |3.5[ it remains to show that 


l(ViV> 


I E l 


ll(Vw) 


(3.41) 


where the mean incorporates the co-normal of the lifted edge E l and the discrete edge E, respec¬ 
tively. Lifting the (Vr h v) to T, equivalence (3.41) can easily proven as in (6.16) (6.18) by using 
the previous lemma to bound the deviation of the lifted discrete co-normal flux from the co-normal 
flux of the lifted edges. □ 


4. Approximation Properties 

Before we turn to the stability and a priori analysis in the next two sections, we here collect 
some standard inequalities and establish the appropriate interpolation estimates which will be 
frequently used throughout the remaining work. 


4.1. Useful Inequalities. First, we recall the following trace inequality for v £ H 1 (Th) 

II^I|9t</i- 1/2 ||h||t + /i 1/2 ||Vi;||t VT sTfc (4.1) 

If the intersection rnT does not coincide with a boundary face of the mesh then the corresponding 
inequality 

IMIrnT</U 1/2 |M|T + /i 1/2 ||Vp|| T VT £% (4.2) 


holds whenever the surface T is reasonably resolved by the background mesh, i.e. the mesh size is 
chosen in accordance with the local curvature, see Hansbo et al. [23] for a proof. Correspondingly, 
since the skeleton T h consists of shape-regular faces, we have 


INI ehf < h- 1 ' 2 |M| f + h l ' 2 \\V7v\\ F VE £ £ h , VF £ T h 
and if we iterate using 4.1 we see that for v £ H 2 {Th) 

IMIsnF < h - 1 |M| t + ||Vp||t + fc||V ® Vti|| r VE £ £ h , VF £ E h 


(4.3) 

(4.4) 


In the following, we will also need the some well-known inverse estimates for Vh £ Vh' 

I|V^||t</i- 1 |Mt VT£Th (4.5) 

Kllar < /U 1/2 K||t, ||Vt> h ||ar < ft" 1/2 ||Vu h || T VT £% (4.6) 

and the following “cut versions” (note that K fl T % dT , E n F dF ): 

ll^lknr < h-^WvhWT, ||V^|knT</i“ 1/2 ||Vu h ||T VK £ K h , VT £% (4.7) 

\\v h \\EnF<h- 1 / 2 \\v h \\ F , \\Vv h \\ EnF < h-^WVvhWr VE£S h ,VF£E h (4.8) 

which are an immediate consequence of similar inverse estimates presented in Hansbo et al. [23] . 
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4.2. Interpolation Operators. Next, we recall from Scott and Zhang [35] that for v £ H m (Th), 
the standard Scott-Zhang interpolant Ih : L 2 (Nh) —> Vh satisfies the local interpolation estimates 


\l-k | 


\\v — IhV\\k,T ^ h "\ v \l,u{T)i 
||w - hv\\i,F < h l ~ k ~ 1 /a \v\i fw(F ), 


0 ^ k ^ l ^ min{2, m} VT S T/i (4-9) 

0O<Z-l/2<min{2,m}-l/2 V F £ F h (4.10) 


where u>(T) denotes the patch of all elements sharing a vertex with T. The patch oj(F) is defined 
analoguously. Before we introduce a suitable interpolation operator Ih : H m (T) -+ Vh, we note 
that by the coarea-formula (cf. Evans and Gariepy [TBI f 


'Us 


f{x) dx = j ^ f(y, r ) dT r (y) j dr 


the extension operator v e dehnes a bounded operator H m (T) 9 v 
the stability estimate 


\ v llfc,i/ 5 (r) 


<5 1/2 I 


't’llfc.r. 


0 sS k € 


v e G H m (Us r)) satisfying 


(4.11) 


for 0 < (5 ^ 5o where the hidden constant depends only on the curvature of T. With the help of the 
extension operator, we construct an interpolation operator Ih : H m (T) —> Vh by setting IhV = IhV e 
where we the liberty of using the same symbol. Choosing 5 o ~ h , it follows directly from combining 
the interpolation estimate (4.10) and the stability estimate (4.11) that the following lemma holds: 


Lemma 4.1. Setting ej = v e — IhV e then for v £ H 2 (T), it holds that 

E ( ft “ 2 |MI F + HVe/ll! + h 2 ||V ® VehWl) < h 2 \\v\\l r 


(4.12) 




The interpolation operator satisfies the following interpolation error estimate with respect to 
the discrete energy norm 


Lemma 4.2. For v £ H 2 (V), it holds that 

|||i’ e - hv e \\\*,h + UK ^ hv e \\\r h < h\\v\\ 2 ,r 

Proof. By dehnition 

IIK - Ihv e \\\l h = (V Fh {v e - I h v e )y Th (v e - I h v e )) Kh 

+ HVr h (v e - I h v e ) ■ n E ,h\\ 2 £h + /i _1 ||v e - hv e \\ £h 

= I+ 11 +III 


(4.13) 


(4.14) 

(4.15) 


which we estimate next. 


Term I. Successively employing the trace inequality (4.1), the interpolation estimate (4.9) and 
the stability bound (4.1 1|) with So ~ h, it follows that 


I= E( V r k (» e -W,Vr k (/^4» e ))r t nT (4.16) 

T&T h 

< E ||VK-4^)11^ (4.17) 

TeT h 

< E (^ _1 l|V(v e - 4i> e )||r + ft||V 8) V (u e - I h v e )f T ) (4.18) 

TeTh 

< E h \\v e \\l MT) (4-19) 

TgTh 

~ ^ l \\ ve \\2,U 6o {T) (4.20) 

<h 2 |M| 2 r (4.21) 
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Term II. To estimate the second contribution, we employ the trace inequaltiy (4.3) to pass from 

(4.22) 


E to F which in combination with (4.10) and (4.11) gives 


■■“'ll 


II < E h\\Vr h (v e - I h v e 

E(z£h 

< E (H V K - IhV e )\\ 2 F + * 2 ||V 8 V(u e - I h v e )\\l) 

E(z.£h 

<h 2 ||«||l,r 

Term III. Applying the same reasoning as for the estimate of Term II, we obtain 

IH< E h ~'\\ e i\\E< E (^ 2 ||e/|| 2 F+ 1^6,111) <h 2 \\v\\l r 

E^Sh E€LJ~h 


To estimate the remaining term |||e/|||.F h , we simply apply the previous Lemma 4.1 
the desired bound. 


(4.23) 

(4.24) 

(4.25) 

to arrive at 

□ 


We conclude this section by briefly reviewing the the Oswald interpolator, cf. Burrnan et al. [4], 
which will be of particular use in establishing discrete Poincare-type estimates in Section (IT2| By 
construction, Oh : “P^JfT) —► V^(Th), with V^ c (Th) and V^(Th) denoting the space of discontinuous 
and continuous piecewise polynonials of order k. Therefore the Oswald interpolator can be thought 
of as a certain “smoothing operator” and the following lemma shows how the fluctuation Id — Oh 
can be controlled in terms of jumps on the faces. 

Lemma 4.3. For v G V^ c {Th), it holds that 

\\Oh{v)\\T < IMU(t) 

\\v-o h {v)f T < E ^IIMIIf 

FeE h (T) 

where Eh(T) denotes the set of all faces F G Th with LOT ^ 0. 

A proof can be found in e.g. Burman et al. [4], together with the construction of the Oswald 
interpolator and a summary of some of its important properties. 


(4.26) 

(4.27) 


5. Stability Estimates 


In this section, we demonstrate that the bilinear form a^(•, •) + jh{-, •) defined in Section 2.4 


is 

both stable and coercive with respect to its associated energy norm. The proof of the coercivity 
result is based on a special discrete Poincare inequality which shows that for discrete functions, 
various L 2 norms can be controlled by merely the tangential gradient and the semi-norm induced 
by stabilization form jh (•,•). As a major instrument in establishing such Poincare inequalities, 
we start by reviewing certain geometric covering relations of the active background mesh and a 
certain stabilization mechanism provided by the bilinear form jh(-, ■). 


5.1. Fat Intersection Covering and Ghost Penalties. Since the surface geometry is repre¬ 
sented on fixed background mesh, the active background mesh Th might contain elements which 
barely intersects the discretized surface r^. Such “small cut elements” typically prohibit the 
application of a whole set of well-known estimates, such as interpolation estimates and inverse in¬ 
equalities, which typically rely on certain scaling properties. As a partial replacement for the lost 
scaling properties we here recall from Burman et al. |6] the concept of fat intersection coverings 

of Th¬ 
in Burman et al. [6] it was proved that the active background mesh fulfills a fat intersection 
property which roughly states that for every element in 7 a there is close-by element which has a 
significant intersection with P/j. More precisely, let i be a point on T and let B$(x) = {y G : 
\x — y\ < <5} and D$ = Bg{x) D T. We define the sets of elements 

ICs, x = {Ke/C h :K l n Ds(x ) ^ 0}, T s , x = {T G % : T n T h £ 1C S , X } 


(5.1) 
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With S ~ h we use the notation and Th,x. For each Th, h £ (0, ho] there is a set of points 
Xh on F such that {ICh, x , x £ X F } and {Th, x -,x £ Xh} are coverings Th and Kh with the following 
properties: 

• The number of set containing a given point y is uniformly bounded 

#{x£X h :y£T h , x }< 1 Vye R d (5.2) 

for all h £ (0, h 0 ) with h 0 small enough. 

• The number of elements in the sets Th, x is uniformly bounded 

#r M < 1 Vx£X h (5.3) 

for all h £ (0, ho] with ho small enough, and each element in T~h, x shares at least one face 
with another element in Th,x- 

• V/i £ (0, ho] with ho small enough, and V x £ Xh, 3 T x £ Th, x that has a large (fat) 
intersection with T^ in the sense that 

\T X \ ~ h\T x nr h \ = h\K x \ Vx£X h (5.4) 

To make use of the fat intersection property, we will need the following lemma which describes 
how to the control of discrete functions on potentially small cut elements can be transferred to their 
close-by neighbors with large intersection by using the face-based stabilization terms appearing in 


Lemma 5.1. Let v be a discontinuous, piecewise linear function defined on a quasi-uniform, mesh 
7ft consisting of element {T} and consider a the macro-element M. = Tj U formed by any two 
elements Xj and T 2 sharing a face F. Then 

IM& < \\vf T2 +h\\[u\f F + h 3 \\n F ■ [Vu]|| 2 f (5.5) 

l|Vv|& < ||Vu||| 2 + /i _1 ||[u]|| F + h\\n F ■ [S/v]f F (5.6) 

with the hidden constant only depending on the quasi-uniformness parameter. 


Proof. A proof of the first estimate can be found in Massing et al. [22]. Setting v = Vw in (5.51, the 
second inequality follows directly from the first one once the following face base inverse estimate 
is established: 


Ml 


< 


| n F ■ [Vw]|| F + h [|[u 


12 
If 


To prove (5.7), rewrite Vu = n F ■ Vu + P F \7v with P F = Id — n F 
face tangential part of the gradient satisfies ||[-PfVu]|| f < h _2 ||[u]|| F . 


(5.7) 

n F and observe that the 

□ 


5.2. Discrete Poincare Estimates. Next, we derive a discrete Poincare estimate showing that 
for v £ Vh the L 2 norm on the active background mesh can be controlled by the tangential 
surface gradient and the semi-norm induced by jh (•,•). While the final Poincare estimate bears 
resemblance to estimates presented in Brenner [T], there are two major differences. First, the 
active background mesh and thus the domain N F varies with the mesh size, and second, the full 
gradient of v £ 14 is not available in our surface method. We start with the following lemma. 


Lemma 5.2. For v £ Vh, the following estimate holds 

h-'Wv-XrMf^ < \\V rh v\\ 2 Th +/i||Vu|| 2 ^ + h~ 2 \\[v]\\ 2 Fh 

for 0 < h < ho with h 0 small enough. 


(5.8) 


Proof. Without loss of generality we can assume that \r h {v) = 0. Apply (5.51 and (4.6) to obtain 

E IMIr hl . (5-9) 


IMIk < 


xex h 

z E 

xCiX h 


Ml k 


< 


Ei 

xGX h 


h\\[v 


2 

Th 


■ h 3 \\n F ■ [Vi>][| F/i 

h 2 \\Vvf Nh 


(5.10) 

(5.11) 
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Thus it is sufficient to estimate the first term in (5.11). For v £ 14, we define a piecewise constant 
version satisfying v\t = pfr f T vdx. Clearly ||u — v\\t < /i||Vu||t- Adding and subtracting v gives 

E imi 

xex h 


E 

<^ 2 ||Vn||^ 


Ml k 


<h 2 \\Vv\\% h 

<^ 2 ||Vn||^ 


E ||w|1 

x&Xh 

E 

x£X h 

^IMIr h + h\\ v ~ ^llr h 

h\\v\\ 2 rh 


(5.12) 

(5.13) 


(5.14) 

(5.15) 


It remains to estimate the last term in (5.15). To apply a Poincare estimate on the surface T, we 
introduce a smoothed version v = Oh(v) £ H 1 (T) of v. Then 

|2 <- J,||~||2 , MU, ~l|2 mi~|| 2 , IU, ~||2 M|~i||2 , mil,1I|2 


h\\vr Th <h\\v 


lr h 


h \\v-v\\ Th < 


Nlr h + ll«-«Hk^^ 


(5.16) 


h \Mw? h 

where we used the inverse estimate (4.7) to pass from T^, to Nh and the fluctation control ( 4.27| ) 
for the Oswald interpolator. Now applying a standard Poincare estimate to the first term yields 


l^'llr < hXriy 1 ) 2 + h\\N r v l \\l =/ + //. 


(5.17) 


which we estimate next. 

Term I. Since Ar h (v) = 0, the first term can be considered as the error of the mean value which 
can be bounded by 

I = h( \ r (v l ) - \ Th (v)) 2 <\v - v\ 2 dT + (1 - cfv 2 dr) (5.18) 

with c = |r /l ||r|” 1 |B|. We note that ||1 — < h 4 thanks to ( |3.23| ). Then 

\v-v\\ 2 Nh +h 4 \\v\\ 2 Nh < hMW^+h^Mk (5.19) 


I < 


h 5 \m\ 2 < 


where the estimates (4.27) and (4.26) for the Oswald interpolator were used in the last step. 
Term II. To estimate the second term, apply the inverse estimates (4.7) and (4.5) to obtain 

TT ^ UIV7_ ^TII2 


1- 1 ||P" 2 


N h 


< 


,<^|Vr^||^+/i||Vr h (u-u)|| 2 h 

(5.20) 

<h\\^r h v\\ 2 rh + \Mv-v)f Nh 

(5.21) 

<mT h v\\l h +h- 2 \\v-v\\l h 

(5.22) 

zH^vwib+h-'mw 2 ^ 

(5.23) 

we see that 


k + HVv\\l h + h-^[v]\\% h +h*\\v\\% h 

(5.24) 


which gives the desired estimates since the last term can be absorbed into the left-hand side for h 
small enough. □ 


The next lemma describes how the full gradient Vu on Nh can be eliminated from (5.8). The 
main idea is to apply the previous lemma to Vu to show that the ft!/ 2 -scaled, full gradient can be 
controlled in terms of the tangential gradient and the face stabilization (2.26). 

Lemma 5.3. For v £ 14, the following estimate holds 


h||Vi 


\N h 


<h 2 \\Vr h v\\k + 


\F h 


(5.25) 


An immediate consequence is 

Corollary 5.1. Let h £ (0,/iq] with ho small enough. Then the following estimate holds: 


h ||u-A r ,» 


\N h 


< 


IIVr.4 


lr„ 


|2 

Wfh 


Mv £ V h . 


(5.26) 
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As a result, ||| • |||^ defines a norm on I4 0 which satisfies the discrete Poincare estimates 


i" 1 lli ; 112 


N h 


< 


< 


where the second version is obtained from the first using the inverse estimate (4.8). 
Proof. (Lemma |5.3| ) Choosing a = Ar h (Vu), the previous lemma gives 

h\\Vv\\ 2 rh < h\\a\\ 2 Th + h\\Vv - a\\ 2 Th 


< 


h\\, 


I T h 


| [Vi 


II Fh. 


It remains to estimate the first term /i||a||j- h . Clearly, 


r h 


h\H 2 Th <h 2 \\a\\ 

and since ||Pra||r fc ^ ll a l|r by a finite dimensionality argument, we obtain 


\ a \\r h < ||«||r < ll-fra||r < ll-Pra||r h ^ ll^r h a||r h + \\{Pr h ~ ^r)a||r h 

^ ll-fVh a llr h + h\\a\\r h 


(5.27) 

(5.28) 


(5.29) 

(5.30) 

(5.31) 

(5.32) 

(5.33) 


We can kick-back /i||a||r fe and absorb it on the left-hand side whenever h ^ ho for some ho small 
enough. Consequently, 

i,2 11 ii ^ j 2 | 


Mr h <h z \\Pr h a\\ 2 rh 


(5.34) 


Now using the boundedness of Pr h and the inverse inequality (4.7), it follows that 
^ll^all^ < ft 2 ||Pr h Vv||p h + h 2 \\P rh (a - Vn)||^ 
<h 2 \\X7 rh v\\ 2 rh + h\\a-Vv\\ 2 Th 


<^ 2 ||Vrw|| 


II Ml 


•Fh 


Collecting (5.31), (5.34) and (5.35) and using ( |5.7[ ), we conclude that 

h\\a\\ 2 n <ti 2 \\Vr h v\\r h 


|| [Vi 


I F h 


<^l|Vr^|lr h + IIK-Vn]||^+fi- 2 ||[u 


| 2 

W^h 


which in combination with (5.30) yields the desired inequality. 


(5.35) 

(5.36) 

(5.37) 

□ 


IMIlk 


(5.38) 


5.3. Coercivity and Continuity. 

Lemma 5.4. The following estimate holds 

h\\n E ■ Vr^lH, < HVr^ll 

for 0 < h < ho with ho small enough. 

Proof. We start with observing that since both nr h and Vw are piecewise constant functions on 
Th so is Vr h u Then successively employing the inverse estimates (4.8),(4.6) 

. |l w ,,||2 


S^IIVr^ 

(5.39) 

<h~' E ||Vr fe n||^ 

(5.40) 

xex h 

<h~' J2 l|Vr^||^ + ||[V r ,i;]||^ 

(5.41) 

xeXh 

< E H v r^||L + ll[Vr,n]|| 2 , ;i 

(5.42) 


x£X h 


where in the last to steps we applied Lemma 5.1 Eq. 5.5 to Vr h and the fat intersection covering 
property as defined in Section|5.l| Since the first term in (5.42) is bounded by ||Vr h w||r , it only 


remains to estimate the second term. Denoting with |p = A( v£ + v T ) the mean value of any 
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piecewise smooth function or operator accross F , standard identities for the jump-operator and 
the boundedness of Pr h give 

ll[Vr^]||^ h < ||[P r J fVv}} \& h + || {{P r J [Vw]||^ < ||[P r J {{Vc}} ||^ + ||[Vu]||^ = I + II 

(5.43) 


To conclude the proof we need to estimate [PrJ- But as in the proof of Lemma 3.1 we see that 

ll[-frJllL°°(F) = \\[Pr h — -Pr]||z,oo( F ) < h 2 (5.44) 

and hence by Lemma |5.3| we arrive at the final estimate 

(5.45) 


Ten. 


□ 


Before we turn to the main stability result, we quickly note as an immediate consequence that 

||MII*,7i ~ HMIU VweVft (5.46) 

Proposition 5.1. The discrete bilinear form Ah{-,-) is both coercive and stable with respect to 
the discrete energy-norm (2.29), that is it satisfies 


\l<A h (v,v) W G 14 


and 


Ah(v, w) < HHIU \\\w\\\ h \/v,w G V h 
whenever /3 E , Pf and 7 are chosen large enough. 


(5.47) 

(5.48) 


Proof. (5.47): Starting from the definition of Ah and applying Cauchy’s inequality with e, we 
have 


A h(v,v)>\\Vr h v\\ 2 K:h -e\\h 1/2 (n E -Vv)\\ 2 £h +(PE-e *))\h %/2 [v]\\ 2 £h + \\\v\\\r h (5.49) 

and employing Lemma |5.4| we immediately see that 

Ah{v,v) > ||Vr h r||| h + \\h~ l/2 [v]\\ 2 £h + |||v||k fe 
when choosing e small and (3 E large enough. 

( 5.48 ): Follows directly from the Cauchy-Schwarz inequality and the norm equivalence (|5.46[). □ 


6 . A Priori Error Estimates 

The goal of this section is to formulate and prove the main a priori estimates for the proposed 
cut discontinuous Galerkin method. We proceed in two steps. First, an abstract Strang-type 
lemma is established which reveals that the overall error can be attributed to two sources, namely 
an interpolation error and a quadrature error caused by the mismatch of the smooth surface and 
its discrete counterpart. Then the quadrature error is bounded using the geometric estimates from 
Section [31 


6.1. Strang’s Lemma. Recalling the definition of the lifted discontinuous Galerkin form (2.34), 
we can state 


Lemma 6.1. With u the solution of (2.7) and Uh the solution of (2.28) it holds 
11 a nJH* /i. "A: |||n Ihu HI* h 


( 6 . 1 ) 


+ sup \Wv\Wb 1 (a h (I h u e ,v) - a l h {I l h u e ,v 1 )) 
vev h v ' 

+ sup IIMIk 1 ^) - l h(vj) 

•uf=LA V / 
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Proof. Let eh = IhU e — Uh • Thanks to triangle inequality |||w e — ^ |||ti e — IhU e \\\* t h~^\w (i h^* 1 h 

and the equivalences of the discrete energy norms (5.46) for v G T4, it is enough to estimate |||e/ lr |||/ 1 . 
We have 


= A-h( K Ih'lL > Cji) Ihiph) 

= A h (I h u e , e h ) - a! h (u,e l h ) + l(e l h ) - l h (e h ) 

= (a h (I h u e ,e h ) - a l h (I l h u e , e l h )^J + ( l(e l h ) - l h (e h ) 

- a l h (u- I l h u e ,e l h ) +j h (lhu e ,e h ) 

= I + II + III + IV 

The bounds for first two terms are immediate: 

I < 111Cfc111h sup IIHII f 1 (a h (I h u e ,v) - a l h (I l h u,v‘)) 
vev h v J 

SU P IlklllhT^) - h{vj) 

vev h v ' 


II < 


Thanks to (3.40) and (5.46), we also have 
III< \\\u - I l h u e \\U 


■'hill* 


< llk e - I h u e 


|*,h 


IV = jh{lhu e — u e ,eh) < |||u 6 - kiW^lkJIkhHk,, 
Collecting the estimates and dividing by |||ek]]h concludes the proof. 


( 6 . 2 ) 

(6.3) 

(6.4) 

(6.5) 

( 6 . 6 ) 

(6.7) 

( 6 . 8 ) 

(6.9) 

( 6 . 10 ) 

□ 


6.2. Quadrature Error Estimates. Next, we show how the abstract quadrature error appearing 
in Strang’s lemma can be bounded in terms of the mesh size and the discrete energy norm. 


Lemma 6.2. The following estimates hold 

\a l h (v l ,w l ) - a h {v,w)\ < ft 2 ||kllkllMIU Vn,w £ V h 
\l(v l ) - l h {v)\ < h 2 \\\v\\\ h \/v £ V h 


( 6 . 11 ) 

( 6 . 12 ) 


Proof. ( | 6 . 1 1 | ) : After lifting the bilinear ah(-,-) from Th to T, each of its contribution can be 
estimated by successively employing the bounds for determinants (3.23)—(3.24), the operator norm 
estimates and the norm equivalences ( |3.25[ )-( |3.26[ ). In doing so we obtain for each K £ Kf 

(ViV, S7tw 1 )k 1 ~ (Vr h v,Vr h w)K 

= (V r v l ,Vw l ) K , - ((V Fh v) 1 , (V^wYlBlf 1 )^ (6.13) 

= ((Pr ~ \B\- L \BB T )W r v\W r w l ) K i (6.14) 

</k||V r k|Uk|ViV|k' (6-15) 


Keeping the estimate for the lifted discrete co-normal (3.30) in mind, each edge contribution gives 
((n E ■ Vrk), [ w l ]) E i - ({n h ,E • Vr h w), M)e 

// w / 


and similarly 


,-lr u 


< 


h 2 \\h 


[w l ]) E i — ((( rih, E • B T X? r v l ), k’1)s 

(6.16) 

1 2 Bn h,E ) • V r k), [w l ]) E i 

(6.17) 

)\\ E ‘\\h~ 1 / 2 [w l ]\\ 2 E i 

(6.18) 

;>lklk" 1/2 k] III 

(6.19) 

'V-(/*">]>])* 

ki-lsi^k^k'D^ 

( 6 . 20 ) 

^ 1/2 k]llllk _1/2 k]lll: 

( 6 . 21 ) 
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(6.12): For the right hand side we have 

Kw l ) - l h (w) = (/, w l ) r - (/ e , w)r h 

<h 2 \\w l \\ r 

< h 2 ||Mlk 


where in last step, the Poincare inequality (5.28) was used after passing from V to T/,. 
6.3. Error Estimates. Finally, we can state and prove the main a priori estimates. 
Theorem 6.1. The following a priori error estimates hold 

UK - Uh\\\*,h < h\\f\\r 
\\n e -u h \\ rh <h 2 \\f\\ r 


( 6 . 22 ) 

(6.23) 

(6.24) 

(6.25) 
□ 


(6.26) 

(6.27) 


Proof. (6.26): Combining Strang’s Lemma 6.1 the interpolation error estimate (4.13), the quad¬ 


rature estimates in Lemma |6.2| and finally the elliptic regularity estimate (2.11) the proof follows 
i mmed iately. 

(6.27): Since the average Ar (u l h ) = |F| 1 f r u l h is not equal to zero we decompose the error as 
follows 


e = u-u l h =u-(u l h - A r (u l h )) + Ar h (u h ) - A r (u l h ) 


(6.28) 


where the first term has average zero on T and the second term is the error in the average. 

To estimate e we let $ £ H 1 ( r)/R be the solution to the dual problem (Vrv, Vr0)r = 
{e,v)r W £ i4(r)/R. Then we have the elliptic regularity estimate ||()>||#- 2 ( r ) < ||e||r• Exploiting 
the fact that [</>] \ l E = [ n E ' 4>\e — 0 by regularity and setting v = e and we obtain 


= ah(e, (f) 

(6.29) 

0) 

(6.30) 

= a l h (e , 0 - ( 40 )*) + a l h (e, (. I h (f>) f 

(6.31) 

< IINII* 1110 ~ (40 e )*|||* + \a l h (e, (40 e )'| 

(6.32) 

< he HI* |0 \\h 2 (t) + a L( e > (40) Z I 

(6.33) 

<h 2 ||/|| r ||e1|r + k(c(4Kl 

(6.34) 


where we added and subtracted an interpolant and estimated the first term using the interpolation 
estimate (4.13) followed by the elliptic regularity estimate for the dual solution. To estimate the 
second term on the right hand side we note that we have the identity 

a(e, v l ) = a(u, v l ) — a(u l h , v) 

= l(v l ) - l h (v) + a h (u h ,v ) - a(u l h ,v l ) \/v £ V h 
Thus using the quadrature estimates (|6.11|)) and (6.12)) we obtain 


| 4 (e,(/K)')l<^ 2 (ll/llr + ||KIIU)|||4r 


(6.35) 

(6.36) 

\hJ\\\*h<P III h < 4||/|| r ||e||r (6.37) 

where at last we used the stability UK 14 < ||/ e ||r h ll/llr of the method and the energy norm 
stability of the interpolant together with energy stability of the solution to the dual problem. 
Together estimates (6.33) and (6.37) give the estimate 

l|e||r<h 2 ||/|| r (6.38) 

To conclude the proof, we note that error in the average e c can be exactly estimated as the last 


term in (5.18), yielding 


I e c11 < |ArK) - ArJu/JI < h 2 |K||r,, < h 2 


(6.39) 
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□ 


7. Condition Number Estimate 

Next, we show that the condition number of the stiffness matrix associated with the bilinear 
form (2.28) can be bounded by 0(h ~ 2 ) independently of the position of the surface T relative to 
the background mesh Th- 

Let be the standard piecewise linear basis functions associated with Th and thus Vh = 

Xu=i Vi4>i for Vh £ Vh and expansion coefficients V = {V.}^L 1 £ R^. The stiffness matrix A is 
given by the relation 

(AV, WOrjv = A h (v h , w h ) V v h , w h £ V h (7.1) 

Recalling the definition of Vh and propositio n|5.1[ the stiffness matrix A clearly is a bijective linear 
mapping A : —> ker(^4)- L where we set R w = M^/ker^) to factor out the one-dimensional 

kernel given by ker.4 = span{(l,..., 1) T }. The operator norm and condition number of the matrix 
A are then defined by 

\\AV\\t 


M||r« = sup ■ 
uei N \o 


N 


and k(A) = ||^1|| r jv||^ 1 || R iv 


(7.2) 


respectively. Following the approach in Ern and Guermond [15] , a bound for the condition number 
can be derived by combining the well-known estimate 

h d/2 \\V\\ RN < |MU 2 dv„) < h d ' 2 ||E|| r « (7.3) 

which holds for any quasi-uniform mesh Th, with the Poincare-type estimate (5.27) and the fol¬ 
lowing inverse estimate: 


Lemma 7.1. Let v £ V^o then the following inverse estimate holds 

IIMIU ^ h- 3/2 \\v\\ Nh 


Proof. First, we note that employing the standard inverse estimates (4.6), we obtain 


M\% h <h- 3 \ 


-Ilk 


(7.4) 

(7.5) 


Now, since Vu| t is constant, an application of the inverse estimates (4.7) and (4.5) gives 


\PvM\l h < ||V«Hk < h- 1 ||v«||k < h~ 3 \\v" 2 


Nh 


(7.6) 


Similarly, recalling (4.8), the co-normal flux and edge-related jump terms can be bounded via 

12 


\\n E ■ Z ^ 2 HV^||k < h~ 3 \ 


\N h 


M\k<h 


- 2 "P " 2 




< h~ 3 I 




(7.7) 

(7.8) 

□ 


which concludes the proof. 

We are now in the position to prove the main result of this section: 

Theorem 7.1. The condition number of the stiffness matrix satisfies the estimate 

k{A) < h~ 2 (7.9) 

where the hidden constant depends only on the quasi-uniformness parameters. 

Proof. We need to bound ||^4||rjv and ||^4“ 1 || r jv. To derive a bound for ||yf|| R jv, we first observe 


that for w £ Vh, 


MIU ^ h- 3/2 \\w\\ Qh < ft. (d - 3)/2 ||W|| R iv 


(7.10) 


where we successively used the inverse estimate ( |7.4[ ) and equivalence ( |7.3[ ). Consequently, 

(7.11) 

Next we turn to the 


M^IIr" = sup = sup ||w|f ^ ~ ^ ^IIMIU ~ ^ ®\ V \ 

II FF || K d wGV h IIIR'IIU IwllRiv 


and thus by the definition of the operator norm, we have ||v4||rn < h d 3 . 
estimate of ||^1” 1 || r jv. Starting from (7.3) and combining the Poincare inequality (5.27) with the 
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stability estimates (5.47) (5.48) and a Cauchy Schwarz inequality, we arrive at the following chain 
of estimates: 

ITII R » < h~ d \\v\\l h < < h l - d A h {v,v) = h}~ d {V,AV) wN < /i 1 " d ||V|| R «||^V|| R A 

(7.12) 

and hence ||V|| r jv < h 1_d || v 4V^|| R N. Now setting V = A~ 1 W we conclude that ||^4 _1 || R JV < 
and combining estimates for ||.4|| r n and ||^4~ 1 || r jv the theorem follows. □ 


8. Numerical Results 

We conclude this paper with two numerical studies. First, we corroborate the theoretical a 
priori estimates presented in Section |6.3| with two convergence experiments. The second study 


serves to illustrate the effect of the different stabilization terms in (2.26) on the sensitivity of the 


condition number with respect to the surface positioning in the background mesh. 

8.1. Convergence Rate. To examine the theoretically expected convergence rates, we consider 
two numerical examples for the Laplace-Beltrami-type problem 


—Arw + u = / on r 


( 8 . 1 ) 


with given analytical reference solution u and surface T = {x £ R 3 : <fi(x) = 0 } defined by a known 
smooth scalar function </> with \7(j)(x) ^ OVa; € T. The corresponding right-hand side / can be 
computed using the following representation of the Laplace-Beltrami operator 


Vu nr — tr(Vnr)Vu • nr 


Aru = An — nr • V i 
For the first test example we chose 

. (irx \ . ( Try \ . / nz \ 

Ui = sm it J sm It) sm It ) 

= x 2 + y 2 + z 2 - 1 
while in the second example, we consider the problem defined by 
U 2 = xy — 5y + z + xz 

fa = ( X 2 - l) 2 + (y 2 - l) 2 + ( Z 2 - if + ( X 2 + y 2 - 4 f + (X 2 + Z 2 - 4 f 

+ ( y 2 + z 2 - 4) 2 - 16 


( 8 . 2 ) 


(8.3) 


(8.4) 


Starting from a structured mesh 7o for fl = [—a, a ] 3 with a large enough such T C fl, a sequence 
of meshes {7fc }| =0 is generated for each test case by successively refining 7o and extracting the 


corresponding active background mesh as defined by (2.14). Based on the manufactured exact 


solutions, the experimental order of convergence (EOC) is then calculated by 


EOC(fc) = 


log(-E fc _i/-Efe) 

log( 2 ) 


where Ek denotes the error of the numerical solution Uk measured in a specified norm and computed 
at refinement level k. In our convergence studies, both || • ||jyi(r^.) and || • ||z, 2 (r fc ) as well as || • ||L=o(./v h ) 
are used to compute Ek. For the two test cases, the resulting errors for the sequence of refined 
meshes are summarized in Table [T] and Table [2j respectively. The observed EOC confirms the 
first-order and second-order convergences rates as predicted by Theorem |6.1| Additionally, we 
observe an optimal, second-order convergence in the L°° (Ab,)-norm. Repeating the convergence 
study with /3e = 0, Pf = 500 and 7 as before yields almost identical error reduction rates for the 
simplified version of (2.25) described in Remark 2.1 (see Table [ 2 ]). The computed solutions are 
visualized in Figure [2] 
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Figure 2. Computed solutions for example 1 (left) and example 2 (right). The 
left part of each solution plot shows the approximation Uh as computed on the 
active mesh Th while the right part displays Uh restricted to the surface mesh ICh- 


Level k 

1 \ u k — **l|jri(r h ) 

EOC 

II u k — u \\L 2 (r h ) 

EOC 

Ifofc - «lft=°(AT h ) 

EOC 

0 

1.27-10° 

- 

2.49-10 -1 

- 

1.44-10 -1 

- 

1 

6.87-10 -1 

0.89 

7.28-10 -2 

1.78 

5.10-10 -2 

1.50 

2 

3.73-10 -1 

0.88 

1.93-10 -2 

1.91 

1.35-10 -2 

1.92 

3 

1.74-10 -1 

1.10 

4.81-10 -3 

2.01 

3.70-10 -3 

1.86 

4 

8.65-10 -2 

1.01 

1.20-10 -3 

2.01 

9.42-10 -4 

1.97 

5 

4.34-10 -2 

0.99 

3.01-10 -4 

1.99 

2.39-10 -4 

1.98 

Table 1. 

Convergence rates for example 1 with Pe 

= Pf = 50 and 7 

= 0.01 


Level k 

ll n fc — u \\H 1 (r h ) 

EOC 

11“*: — n lli 2 (r h ) 

EOC 

ll“fc — u \\L°°(N h ) 

EOC 

0 

2.21-10 1 

- 

1.48-10 1 

- 

4.56-10° 

- 

1 

7.33-10° 

1.59 

1.84-10° 

3.01 

8.54-10 -1 

2.42 

2 

3.12-10° 

1.23 

5.41-10 -1 

1.77 

2.26-10 -1 

1.92 

3 

1.61-10° 

0.95 

1.16-10 -1 

2.22 

5.38-10 -2 

2.07 

4 

7.65-10 -1 

1.07 

2.86-10 -2 

2.02 

1.28-10 -2 

2.07 

5 

3.78-10 -1 

1.02 

7.18-10 -3 

2.00 

3.32-10 -3 

1.95 


Level k 

ll“fc — “IlnUr^) 

EOC 

11“* — “lli 2 (r h ) 

EOC 

ll“fc - “lft=°(AT h ) 

EOC 

0 

2.20-10 1 

- 

1.49-10 1 

- 

4.59-10° 

- 

1 

7.22-10° 

1.61 

1.80-10° 

3.05 

8.20-10 -1 

2.48 

2 

3.09-10° 

1.22 

5.38-10 -1 

1.74 

2.23-10 -1 

1.88 

3 

1.62-10° 

0.93 

1.16-10 -1 

2.22 

5.27-10 -2 

2.08 

4 

7.64-10 -1 

1.09 

2.85-10 -2 

2.02 

1.25-10 -2 

2.07 

5 

3.78-10 -1 

1.02 

7.14-10 -3 

2.00 

3.29-10 -3 

1.93 


Table 2. Convergence rates for example 2 with (3e = Pf = 50 (top) and Pe = 0, 
Pf = 500 (bottom) and 7 = 0.01 in both cases. 


8.2. Condition Number Tests. Finally, we numerically examine the mesh-size dependency of 
the condition number of our proposed method and study how the positioning of the surface in the 
background mesh affects the condition number. 















22 A CUT DISCONTINUOUS GALERKIN METHOD FOR THE LAPLACE-BELTRAMI OPERATOR 


Let {7X-}fe=o a sequence of successively refined tessellations of Q = [—1.6,1.6] 3 with mesh 
size h = 3.2/5 • 2~ k . On each refinement level k, we generate a family of surfaces 
by translating the unit-sphere S 2 = {a; £ R 3 : ||x|| = 1} along the diagonal (h,h,h); that is, 
r a = S 2 + So(h, h, h) with <5 € [0,1]. For 6 = 1/ 500, l = 0, ...,500, we compute the condition 
number Kg(A) as the ratio of the absolute value of the largest (in modulus) and smallest (in 
modulus), non-zero eigenvalue. For each refinement level k, the resulting condition numbers are 
plotted in Figure[3]as a function of S. We observe that the position of T relative to the background 
mesh 7 ~k has very little effect on the condition number when the stabilization parameters and 
7 are chosen sufficiently large. In contrast, the condition number is highly sensitive and clearly 
unbounded as a function of 5 if we set either penalty parameter in (2.26) to 0 as the corresponding 
plots in Figure[3]reveal. Finally, scaling the largest computed condition number on each refinement 
level k with k~ 2 confirms the theoretically proven 0(h~ 2 ) bound, see Table [3j 

To further elucidate the importance of the various stabilization terms in ( |2.26 ), we conduct 
a second numerical experiment inspired by the results presented in Olshanskii and Reusken [3T| . 
Olshanskii and Reusken m show that diagonal preconditioning yields robust condition number 
bounds and cures the discrete system from being severely ill-conditioned when a continuous piece- 
wise linear ansatz space is used, see also Figure [3] Thus, no stabilization is needed and 7, as the 
only relevant penalty parameter in the continuous case, can be set to 0. Repeating the same exper¬ 
iment for our proposed cut discontinuous Galerkin method with either penalty parameter /3p or 7 
deactivated shows that the same conclusion is not true if discontinuous finite element functions are 
employed. The reason lies in the construction of the approximation space by restricting the finite 
element space defined on the background mesh to the surface. For each element T £ Th, restricting 
the shape functions defined on T to the surface part K = fl T clearly yields a set of locally 
linearly dependent functions on K. With the basis functions stretching over several elements, this 
effect is usually counteracted in the case of continuous finite element functions, since curvature 
effects lead to different linear dependencies on each element. Thus the resulting set of finite ele¬ 
ment functions on ICh is usually only “nearly” linearly dependent and therefore a well-conditioned 
system can be obtained by means of proper preconditioning. However, in the discontinuous case, 
restricting the finite element basis on 7/ to the surface produces a globally linearly dependent set 
of discrete functions, even in the presence of highly curved surfaces, as interelement continuity is 
imposed only weakly via (2.26). We also observe that if we enforce nearly inter-element continuity 
of the discrete functions by choosing a very large penalty parameter, e.g., /3f = 5e6, diagonal 
preconditioning yields robust, albeit large, bounds for the condition number, see Figure [3] 


k 12 3 4 

r 2 maxj{Kj(i)} 7370.36 7167.30 7205.85 6271.44 


Table 3. Maximum of scaled condition numbers for each refinement level k. 
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Figure 3. Condition numbers plotted as a function of the position parameter 5. 
Top: With either (3p = 0 or 7 = 0, the condition number is highly sensitive to 
the surface positioning in the background mesh. With both penalties activated, 
the condition number is robust and the edge-based penalty term /3e can be omit¬ 
ted. Bottom: Diagonally preconditioned surface CutFEMs. For an unstabilized 
continuous Galerkin (cG) method, preconditioning gives robust condition num¬ 
bers while for our discontinuous CutFEM (dG), the preconditioned system is still 
highly dependent on the surface position if either /3p = 0 and 7 = 0 . 
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